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Executive Summary 

An effective and robust interface element technology able to connect independently 
modeled finite element subdomains has been developed. This method is based on the use of 
penalty constraints and allows coupling of finite element models whose nodes do not coincide 
along their common interface. Additionally, the present formulation leads to a computational 
approach that is very efficient and completely compatible with existing commercial software. A 
significant effort has been directed toward identifying those model characteristics (element 
geometric properties, material properties and loads) that most strongly affect the required penalty 
parameter, and subsequently to developing simple “formulae” for automatically calculating the 
proper penalty parameter for each interface constraint. This task is especially critical in 
composite materials and structures, where adjacent sub-regions may be composed of 
significantly different materials or laminates. This approach has been validated by investigating a 
variety of two-dimensional problems, including composite laminates. 



1. INTRODUCTION 


The ways in which analysis and design are performed have changed extensively during 
the past decade. Automated design algorithms are now plentiful and increasingly robust, and no 
longer are individual components of a structure designed in a vacuum. Facilitated in part by 
product data management (PDM) and product lifecycle management (PLM) software and 
internet-based data sharing, integrated design activities are now possible. Further, the speed and 
economy of modern computers have enabled engineers to perform many large scale analyses that 
were unheard of a decade ago. It is now common to analyze and simulate the response of entire 
aircraft, spacecraft, automobiles, ships, and other structural assemblies to a variety of complex 
and combined types of loading. The CPU time required to perform such analyses is very small 
compared to the time required for engineers to create the mathematical and computer models, 
and the latter effort is usually the most expensive component of a large scale analysis or 
simulation. 

With model sharing and large scale analysis activities on the rise, it is becoming evident 
that improved technology for building computer models is needed. One issue that arises often is 
the need to perform a unified analysis of a structural assembly using sub-structural models 
created independently. These sub-structural models are frequently created by different engineers 
using different software and in different geographical locations, with little or no communication 
between the teams of engineers creating the models. As a result, these models are likely to be 
incompatible at their interfaces, making it very difficult to combine them for a unified analysis of 
the entire assembly. 

Finite element interface technology has been developed during the past decade to 
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facilitate the joining of independently modeled substructures. Unconventional approaches have 
been employed to connect special elements based on analytical solutions to finite element models 
[1-2]. In [1] a near-field solution for a dynamically propagating crack at the interface of two 
dissimilar anisotropic elastic materials was successfully implemented into a hybrid-displacement 
finite element formulation. In [2] Lagrange multiplier terms are used to couple a special element, 
based on a stress analytical solution, to standard finite elements. 

In order to take advantage of parallel computing, Farhat and colleagues [3-4] developed a 
domain decomposition approach for partitioning the spatial domain into a set of disconnected 
subdomains, each assigned to an individual processor. Lagrange multipliers were introduced to 
enforce compatibility at the interface nodes. In [5] non-conforming “mortar” elements are 
employed to connect incompatible subdomains using a conjugate gradient iterative technique in a 
domain decomposition scheme designed for parallel computers. 

The finite element interface technology developed at NASA LaRC [6-10] and elsewhere 
[11] allows the connection of independently modeled substructures with incompatible 
discretization along the common boundary. This approach has matured to a point that it is now 
very effective. However, because the interface technology utilizes Lagrange multipliers to 
enforce the interface constraint conditions, the resulting system of equations is not positive- 
definite. Hence, state-of-the-art sparse solver technology cannot be utilized with the interface 
technology. 

Recently, an alternative formulation for the finite element interface technology based on 

Lagrange multipliers has been developed [12]. The alternative approach recasts the interface 

element constraint equations in the form of multi-point constraints. This change allows an easier 

implementation of the formulation in a standard finite element code and alleviates the issues 
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related to the resulting non-positive definite system of equations. The method seems to provide 
reliable results, but the formulation of the interface method is still quite complicated. 

A possible remedy for these shortcomings is to modify the current hybrid variational 
formulation of the interface element by enforcing the interface constraints via a penalty method 
as opposed to the current Lagrange multiplier approach. The primary consequences of this 
modification will be (i) a simple formulation that is easily implemented in commercial finite 
elements codes, (ii) a positive-definite and banded stiffness matrix and (iii) a reduced number of 
DOFs, since the Lagrange multiplier degrees of freedom (DOFs) will not be present. Thus, the 
proposed approach should greatly enhance the computational efficiency of the interface element 
technology. 

From an accuracy point of view, the penalty method enforces the constraints only 
approximately, depending on the value of the penalty parameter chosen, while the Lagrange 
multiplier approach enforces the constraints exactly. The penalty method interface approach was 
recently attempted using a single global value of the penalty parameter to enforce all constraints 
[13]. This study demonstrated the validity and the effectiveness of the penalty approach in an 
interface element. However, there is need for specific guidelines regarding the selection of an 
appropriate value of the penalty parameter, especially when the substructures to be connected 
have different material and/or section stiffnesses. 

There is a large body of literature related to the determination of an optimal value of the 
penalty parameter (see, e.g. [14-26]). However, there is no existing criterion for choosing the 
penalty parameter in the framework of the interface element under investigation. The 
determination of such a criterion is the primary thrust of the effort presented herein. 

An effective and robust interface element technology has been developed using the 
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penalty method. A significant part of the effort has been directed toward identifying those model 
characteristics that most strongly affect the required penalty parameter, and subsequently to 
developing simple “formulae” for automatically calculating the proper penalty parameter for 
each interface constraint. The new approach has been validated through a wide variety of one and 
two-dimensional problems that have been investigated extensively from both an analytical and a 
computational point of view. Finally, the penalty based interface element technology has been 
implemented into an existing commercial code. 
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2. GENERAL DESCRIPTION OF THE INTERFACE ELEMENT 


Consider two independently modeled subdomains Q, and Q 2 as shown in Figure 1(a) 
and 1(b), respectively, for a 2D and for a 3D geometry. The two substructures are connected to 
each other using an interface element acting like “glue” at the common interface. The interface 
element is discretized with a set of nodes that are independent of the nodes at the interface in 
subdomains Q, and Q 2 . Both in the hybrid interface method and in the penalty method, the 
coupling terms associated to the interface element are arranged in the form of a “stiffness” matrix 
and assembled with the other finite element stiffness matrices as usual. 

2. 1 Hybrid Interface Method 

The hybrid interface method [6-12] introduces two vectors of Lagrange multipliers /, 
and /L, in the total potential energy (TPE) of the system to satisfy the displacement continuity 
conditions. Thus the TPE of the system assumes the form: 

7 Z = x ni +x ai + fA, (V -u i )ds+ -u 2 )ds ( 1 ) 

S s 

The nodal displacements of the sub-domain Q, are identified by q" and q' r The 
superscript o identifies the degrees of freedom (DOFs) that are not on the interfaces, while i 
denotes DOFs that are on the interfaces. The displacement field u t of the sub-domain Q f is 

expressed in terms of the unknown nodal displacements q) as: 

u j= N j<Ij (2) 

where N; can be the matrices of linear Lagrange interpolation functions. 

The displacement field V is approximated on the entire interface surface in terms of 
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unknown nodal displacements q s as: 


V = Tq s (3) 

where T is a matrix of cubic spline interpolation functions [27-29]. 

The first variation of ;ris taken with respect to all the DOFs and the vectors of Lagrange 


multipliers A ] and A 2 : 

8 k L , „ , . . = 0 

On the interface part of the subdomains the following equations result: 

Uj = V, A j = t j} ^+^=0 for 7 = 1,2 

where tj is the traction on the interface. These equations show that: 

• Displacement continuity is enforced, 

• A f are interface tractions, 

• the total traction is zero on the interface. 

The first variation of ^-yields the system of equations: 


(4) 

(5) 


' K" K[° 0 0 0 M, 0 1 f q[ 1 \ff 

K"' K"° 0 0 0 0 0 q\’ f" 

0 0 K“ K'° 0 0 M 2 q!, fi 

o o k; 1 k; 0 o o o <q'; >=</“> (6) 

0 0 0 0 0 G, G 2 q s 0 

Ml 0 0 0 G, r 0 0 \ 0 

0 0 Ml 0 G\ 0 oJ[^j [o 

It can be seen that the resulting global stiffness matrix is sparse, symmetric, not banded and not 

positive definite. Inside [A!] the interface element is represented by the coupling terms Mj and Gj 

which augment the stiffness matrix of the subdomains. 

A “stiffness” matrix and generalized vector of unknown displacements can be associated 


with the interface element, so that: 
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( 7 ) 


2. 2 Penalty Hybrid Interface Method 


In the penalty hybrid interface method the displacement continuity constraint is imposed 
in a least squares sense through two vectors of penalty parameters y and y 2 . Thus the TPE of 
the system assumes the form: 


* = *o, + *h 2 + tTi \{V ~ m, f ds + ~y 2 \{V - u 2 f ds 


( 8 ) 


The displacement fields u l and V are approximated in a similar way as in the original hybrid 


interface method. 


The first variation of n is taken with respect to all the DOFs, but not the vectors of 


penalty parameters y l and y 2 , which are predetermined constants. 

5k L , „ , =0 


( 9 ) 


We now define: 


1 

IK n M 

( 10 ) 



o“ = r , J 

KT)* 

( 11 ) 



II 

A 

(12) 

G J = rj : 

\{ T J T ,) ds 

( 13 ) 


s 
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So the global system of equations of the penalty hybrid interface method assumes the following 
form: 

~K™ Kf 0 o olU°] f/f ' 

k;° k;'+g f -Gf o o q \ /; 

0 -Gf Gf + Gf -Gf 0 < q s l = < 0 l (14) 

0 0 -Gf Kf+Gf K; q[ /' 

0 0 0 Kf K"”\[q"\ [f;\ 

This is a symmetric, banded and positive definite (after boundary conditions are imposed) global 

stiffness matrix. The “stiffness” matrix and generalized vector of unknown displacements 

associated with the interface element can be defined such that: 

" Gf -Gf o ]U;| [o' 

-Gf Gf + Gf -Gf «,Uo (15) 

0 -Gf GfJ^J [o 
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3. DETERMINATION OF THE PENALTY PARAMETERS 

In the penalty method, the displacement continuity constraint is imposed through penalty 
parameters, a set of predetermined constants. The FE solution obtained using this method is 
approximate, with its accuracy depending on the value of the adopted penalty parameters. It is 
known that the penalty parameter should depend on the material and/or geometric properties of 
the two sub-regions being joined. Further, there is a relationship between the penalty parameter 
and the Lagrange multiplier that enforces a given constraint. 

The Lagrange multiplier method imposes the continuity constraint exactly. Thus, the 
Lagrange multiplier method defines the upper limit to the accuracy of the penalty method. 
Knowledge of the correct solution facilitates relating the value of the penalty parameter to the 
geometrical and material properties of the model under consideration. In our pursuit of the proper 
penalty parameter values, a variety of one and two-dimensional problems have been studied with 
both the Lagrange multiplier method and the penalty method. 

The types of finite elements that have been investigated are: conventionally formulated 
and reduced integrated Timoshenko beam elements, plane stress quadrilateral elements and plate 
elements based on the first order shear deformation theory (FSDT), or Mindlin plate theory. For 
each finite element formulation, different penalty parameters are associated to the various nodal 
DOFs. For example, the Timoshenko beam element has three independent nodal DOFs: the axial 
displacement u, the transverse displacement w and the rotation y/. Thus, three different penalty 
parameters y u , y w and y are employed to enforce the interface continuity constraints on the 

DOFs u, w and y/. An independent choice of the penalty parameters is of fundamental importance 
since each degree of freedom can be related differently to the material and geometric properties 
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of the finite element model. 


The methodology adopted in finding the relations will now be described. First the most 
common load cases for the FE type under consideration are applied separately to a simple model 
of one or two elements. For example, for the Timoshenko beam element the following load cases 
are considered: axial load, transversely distributed load, concentrated transverse load and 
concentrated moment applied at the tip. These loads are applied to two different models: a single 
beam element connected to a fixed point by one interface element (see Figure 2a), and two linear 
beam elements connected using an interface element and clamped at the tip (see Figure 2b). 

The formulations and solutions are obtained using both the Lagrange multiplier method 
and the penalty method. The displacement solutions from the two methods are compared 
individually for each degree of freedom. The ratio between the two solutions is expressed in the 
form: 


penally 
y Lagrange 



(16) 


where /=/ (element geometric properties, material properties, and loads) 

Once this simple expression has been identified, the penalty parameter y is set equal to: 

r = Pf ( 17 ) 

Then, the ratio between the solutions becomes independent of material and geometrical 
properties of the element: 


y penally | 

— / — 1 H 

ylMgrange p 


(18) 


The accuracy of the solution depends directly on the value assigned to the parameter f3 . For 
example, if /? is equal to 1000, the penalty solution differs from that obtained by use of the 
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Lagrange multipliers by 0.1%. The degree of precision of the solution cannot be indefinitely 
increased, since round off amplification error would rise. However, once a reasonable 
compromise between constraint representation error and the round off error has been evaluated, a 
value of (3 can be identified that is able to produce the same level of accuracy for every 
combination of material and geometrical properties. 

In some cases, however, the displacement solution related to a particular DOF depends on 
the penalty parameter used to enforce a different DOF. An example of this situation is 
represented by the following relations: 


penally 

u \ 

La grange 

U\ 


= 1 +— +— 
h r 2 


(19) 


u Pwa!ty 

ylMgrange 



( 20 ) 


where m, and u 2 are displacement solutions computed at the same node for two different DOFs, 
while /, , f 2 and / 3 are independent functions of the model properties and the loads. 

Two possible options are available in this situation. The first option is to choose 
Y\=ipf s and y 2 = 2 /?/ 2 in order to obtain the desired simple relation between the solution ratio 


and (3 only for the DOF 1 . 


u penalty 


Lagrange 

U \ 


i 1 1 i 1 

1 H 1 — 1 H — 

2(3 2/3 (3 


y penalty 

Lagrange 
U 2 



(21) 


The other option is to set y , =2 and y 2 = fif 3 , favoring a simple relation for the DOF 2. 


u 


penalty 
1 


iMg range 


= 1 + 
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■ + 


/ 2 


penalty 
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PI: \ 
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1 + I 

Lagrange p 


( 22 ) 


Clearly a conflict arises about the value of the penalty parameter y 2 . A reasonable way to 
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address this issue is to compare the two possible values for y 2 and to choose the greater. Thus, 
the accuracy for one DOF is directly proportional to /? as desired, and the precision of the other 
one should be slightly higher. 

As discussed previously, one of the main objectives of this study is to establish an 
automatic choice of the optimal penalty parameter. The investigations of the ratio between the 
solutions based on Lagrange multipliers and the penalty method constitute the basis for pursuing 
this goal. Nevertheless, it should be underlined that an exact value of the penalty parameter is not 
required. Rather, a value that is of the right order of magnitude is sufficient. In fact, even in the 
most complex FE analysis, there exists a range of values for this parameter for which the 
numerical outcomes change very little. This range can equal as much as 12 orders of magnitude 
for simple analyses, but usually is not less than two orders of magnitude in most situations. 

3. 1 Computational Example -Beam Under Uniaxial Load 

In order to limit the length and the complexity of the mathematical derivation, the 
example described here is one of the simplest available: the extension of a beam under a uniform 
load P applied at the tip. This case is analyzed using a single linear beam element connected to 
one fixed point V by a displacement continuity constraint, which is imposed through a Lagrange 
multiplier or a penalty parameter. The configuration of the geometry and the mesh is plotted in 
Figure 2a. 

3.1.1 Lagrange Multiplier Method 

The hybrid interface method introduces a Lagrange multiplier A t in the total potential 
energy (TPE) of the system to satisfy the displacement continuity condition. Thus the TPE of the 
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system in this study assumes the following form. 


n 


= n x +A, • (v - w, ) - r, -v 


( 23 ) 


where /r, is the TPE of the bar; r t and v are, respectively, the reaction force and the 
displacement at point V. Expanding and simplifying: 

n= | Ua x -£ x )d V -Yjf>' u i + ^\ v -U')- r x- v ( 24 ) 


/=t 


7 T = ~ \ea(^- ] dx-^f, ■u l +A i (v-u ] )-r i -v 
2 n J \dx) , =1 


( 25 ) 


where E is the Young’s modulus, A is the cross sectional area and fi are the applied nodal forces. 
If the displacements are approximated linearly by: 


u = Y J U , N J = M i 
,/=i 


HMf 


( 26 ) 


n takes the following form. 


EA 'WV dN / \ , v 


o \ ' 


k J 


( 27 ) 


Setting to zero the first variation of n , taken with respect all the DOFs, produces the following 


equations. 


dn 

du. 


= EA- J 


L/ dN^ r 


dx 


V dN > 

vT J * 


a* _ \(dN^ r 


dn. 


■=EA-\\ 


dx 


dN. 


y«p 

KJ J J* 


dx -A, = 0 


dx- P = 0 


dTT . „ 

— = A, - r, = 0 
dv 


( 28 ) 


( 29 ) 


( 30 ) 
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( 31 ) 


v - u, =0 


Thus, the resulting FE model is: 


EA 

EA 

L 

L 

EA 

EA 

L 

L 

0 

0 

-1 

0 


0 -1 

0 0 


Since V is fixed, v = 0 , and the solution is found to be: 

u, = v = 0 


The solution provided by the hybrid interface method coincides with the exact theoretical one. 

3.1.2 Penalty Method 

In the penalty method the displacement continuity constraint is imposed through a 
penalty parameter y ] . Therefore the TPE of the system takes the form: 


* = *i + rK '( V-M i) _i 


Tt^lUcr^e^dV-^fru^-y, {v-urf -r t 


1 f r- 4 f du I J r 1 / \2 

= - J EA — dx-Xfi-Ui+Th -( v -»i) -^1 • 


Approximating the displacements linearly: 


15 



( 39 ) 


EA rf v dN.V 

* = — ■ J 2> 


2 i nr ' dx 


^ dNA 1 \ 2 

> w, — - dx-P u 2 +-/] '(v-m,) -r, -v 

V Y ^ J 2 v ' 


The first variation of n in this case is taken with respect all the DOFs, but not the penalty 
parameter. 


dn „ . "ri dN , 
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(41) 

r i = o 

(42) 


The resulting FE model is: 
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Solving the system of equations: 


r] =-P 

P 

M, = — 

Y\ 


r ~ + — 
Yi + EA 


(43) 


(44) 

(45) 

(46) 


3.1.3 Comparison Between the Two Methods 

When the penalty parameter approaches infinity, the penalty method solution tends to the 

one obtained using the Lagrange multiplier method. 
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( 47 ) 


P 

Limit w, = Limit — = 0 
n-+® ri->® y , 


Limit u 2 = Limit 


( 1 

L } 


L } 

— + — 

P = \ 

— 

u 

EA J 


,ea) 


(48) 


Moreover, the constraint that has been enforced is: 

C i (w,,m 2 ,v) = (v-m i ) (49) 

If we define u ir ,u 2y ,v r as the solutions derived from the penalty formulation, it is possible to 
verify the well-known relation between the Lagrange multiplier and the penalty parameter. 

P 


K = 7i •Q( H ir’*V v r) = ?v( V V _ *v) = ^, 


0 — 


7\ 


= -p 


(50) 


3.1.4 Relation Between Penalty Parameter and Beam Properties 

The exact displacement of the tip of the beam matches that from the Lagrange multiplier 

finite element formulation. 


u 


exact 

2 



(51) 


The Penalty parameter solution u$ e " ally differs from the exact one by the presence of an additional 
term (P / y t ): 


The ratio 
solutions: 
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penalty 


_L -L 

U + EA 


(52) 


^penalty j ^exact cva ] ua f ct ] jn order to identify the relationship between the two 
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The penalty parameter y ] is now substituted by: 


(53) 


Y\ =P' 


EA 


(54) 


It follows that the ratio between the solutions becomes independent of material and geometrical 


properties of the beam. 


^ penalty 



(55) 


A similar procedure is followed for the two-dimensional and three dimensional problems 
of any kind. Details are omitted for brevity. But the final results of this investigation are 
presented in the appendix for selected element types. 

A common approach is to set the penalty parameter equal to the largest diagonal term in 
the stiffness matrix multiplied by a factor 10 w , where n is a preselected integer value. For the 
simple example above, Eq. (54) yields an equivalent result. However, in general the approach 
developed here in yields expressions that are not directly proportional to coefficients in the 
stiffness matrix. 


3,2 Automatic Round-Off Error Control 

Due to finite precision in floating-point arithmetic used when the interface element 
stiffness matrix is numerically integrated, the stiffness coefficients are always approximated. 
However, in order to be imposed correctly (and to contribute no energy to the system), the 
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displacement continuity constraint {V — u) requires the sum of the terms in every row of its 
stiffness matrix (15) to be zero. This condition usually cannot be achieved, due to the round-off 
error, and the resulting inaccuracy grows with the value of the penalty parameter. Precisely, the 
important measure is the ratio between the order of magnitude of the interface element stiffness 
matrix rows’ imbalance and the element stiffness. If K n is the stiffness associated to the /7-th 
nodal DOF, it is sufficient to consider the ratio: 


Q=^ 

" K., 


(56) 


where ER is the unbalance in the interface eleme nt stiffness matrix row related to the DOF n. 


ER„ = 'Z k -i (57) 

./ 

When the value of Q„ exceeds about MO" 4 , errors in the solution may become 
appreciable. The discussed row imbalance is proportional to the value ol the penalty 
parameter, ER n oc y . It is also approximately true that: 

ER n oc y x /3 ■ K n => ( 58 ) 

Accordingly, an algorithm has been developed to control the round-off error. Its steps can 
be summarized as follows: 

• Stiffness terms for every nodal DOF in the interface element are computed from known 
geometrical and material properties. 

• For each row in the stiffness matrix: 

o The highest stiffness term is selected and assigned to a variable K 
o The row imbalance of the stiffness matrix is stored in a variable ER 
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o O = is evaluated 

K 

• The highest Q found is compared to a given constant value C . Typically C = 1 • 1 O' 7 is used. 

C 

• If Q > C , the parameter (5 is reduced according to: J3 new = — ■/? 

• The interface element stiffness matrix is recalculated using the new value (3 = f3 m " . 

This approach reduces the risk that round-off errors could adversely affect the solution. 
Thus, the initial value of (3 can be increased, in order to get a higher degree of accuracy, 
knowing that it will be automatically reduced if rounding errors don’t allow that precision to be 
realized. 


3.3 Implementation of the Model as an Abaqus User Element Subroutine 

To test the behavior of the penalty method based “interface” element and the accuracy of 
the obtained relations for an automatic choice of the penalty parameters, the element has been 
implemented in the commercial finite element code ABAQUS as a User Element Subroutine 
(UEL) [30], 

The UEL subroutine receives all the necessary information about geometry and material 
properties of the two connected meshes of finite elements from the input file. Then, the stiffness 
matrix of the interface element is built as previously defined: 


G," 


-G? 


-G," G," + Gf 


0 


-G" 


-G 2 " 

G" 


(60) 


The appropriate values of penalty parameters for each constraint are assigned automatically 
inside the subroutine. 
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4. NUMERICAL RESULTS 


4. 1 Isotropic 2D Cantilever Beam with a Vertical Interface 

A two dimensional cantilever beam is assumed to be composed of two domains. The 
domains are meshed independently and joined by one interface element. The geometrical and 
load configuration of the problem are represented in Figure 3. Others properties are: P = 1, 
Q = 1 , £, = E 2 = 1 and thickness = 1 . 

Only one interface element is adopted in this simple case, since more would not improve 
the solution accuracy. This element is identified by the vertical interface that contains the three 
DOFs V„ V 2 and V 3 . In discretizing the two subdomains, four-noded bilinear plane stress 
elements are used. Two different load conditions are investigated: axial and bending loads. This 
problem may be considered a patch test for the element developed herein. 

4.1.1 Axial Load 

Under the condition of a uniformly distributed axial load applied at the free end of the 
beam, the interface element results are in agreement with the exact solution to the number of 
significant digits available. Results from two different beam discretizations are reported in Table 
1. 

It is important to study the stability of the solution obtained by the new method for 
various values of the nondimensionalized penalty parameter y / E . For this test the automatic 
choice of the optimal parameter y was temporarily disabled. Results are produced in graphical 
form in Figure 4. It can be observed that the solution is stable for a wide range of values of the 
nondimensionalized penalty parameter y/E. 
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4.1.2 Bending Load 

For the bending load case, a comparison is made to both the classical beam solution and a 
reference finite element solution. The coarse mesh of bilinear quadrilateral elements cannot 
exactly recover the classical solution. In Table 2 are results obtained from two different beam 
discretizations. The FE results, marked as Abaqus, are obtained from a traditional analysis using 
a compatible finite element model of the beam. 

The stability of the solution obtained is tested also in this load case for various values of 
the nondimensionalized penalty parameter y/E. Again the automatic choice of the optimal 
parameter y was temporarily disabled. Results are produced in graphical form in Figure 5. The 
method proves very reliable and stable for bending loads. 

4.2 Tension-Loaded Plate with a Central Circular Hole 

A plate with a central circular hole is subjected to a uniform load at its edges. This well 
studied elasticity problem allows one to investigate the capability of the interface elements in 
performing global/local analysis. Geometrical and load configurations are plotted in Figure 6. 
The plate height and width are, respectively, 6 = 100 and 6 = 200; the radius of the hole is 
r = 2.5 ; the distance of the interface from the center of the hole is 5. 

Taking advantage of symmetry, only one quarter of the plate is modeled. A fine mesh is 
applied in the area between the hole and the interface, while a much coarser mesh is adopted 
elsewhere. The complete finite element model is reported in Figure 7a. For a better view, the area 
near the hole is depicted in Figure 8a. The elements used to mesh the two domains are plane 
stress bilinear quadrilaterals, Abaqus CPS4. Five interface elements are employed along each 
segment of the interface for a total of 1 1 interface nodes. 
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Figures 9a and 10a demonstrate that only light discontinuities in the values of stresses at 
the interface are present. Horizontal and vertical displacements do not show any change in their 
values across the interface, as plotted in Figures 1 la and 12a. 

In order to compare the previous maps of stresses and displacements, a conventional FE 
model has been built with Abaqus. Figures 7b and 8b show the adopted mesh. Stresses <r n , c 22 
and displacements U l and U 2 are provided in Figures 9b to 12b. Comparing the maps of the 
displacements, it can be noticed that the interface elements do not introduce any discontinuities 
at the interface. The distribution of the stresses shows small differences, but they can still be 
considered sufficiently accurate. 

The elasticity solution for an infinite plate with a central hole loaded in uniform tension 
predicts that the stress concentration factor, K t , is equal to 3.0 at the upper edge of the hole. 
Moreover, it provides the equation for evaluating the change in a u along the line axis X 2 . 
Using this expression, we computed the stress variation along X 2 for our geometrical 
configuration and compared this result with our FE models. Figure 13 graphs the exact elasticity 
solution against the penalty interface and the conventional FE ones. The results are all very close, 
validating the precision of the developed method. 

4.3 Composite 2D Cantilever Beam with a Vertical Interface 

A 2D composite cantilever beam made of five materials is considered. The geometrical 
and loading configurations are the same as for the cantilever beam analyzed previously (see 
Figure 14). The entire left domain of the beam is made of one isotropic material having a 
Young‘s modulus equal to 2.9E+7 MPa and Poisson’s ratio 0.25. The right domain is composed 
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of four isotropic layers of equal thickness. They have all the same Poisson ratio 0.25, but the 
Young’s modulus E varies. Starting from the bottom, the layerwise E takes the values: 1E+6 
MPa, 1E+7 MPa, 5E+6 MPa and 1E+8 MPa. A conventional finite element model (10x4-10x4) 
has been used as a baseline to test the results from the penalty interface method. 

Under these conditions the interface is expected to undergo abrupt changes of slope. An 
important property of the interface elements is the automatic choice of the penalty parameters. In 
fact, the interface elements can join the four different layers to the left domain with the dissimilar 
stiffness. This behavior can increase the accuracy of the results, avoiding the possible corrupting 
effect of an overestimated penalty parameter. 

The results are shown for two possible interface discretizations: 1 interface element and 4 
interface elements. Axial and vertical displacements at the free end of the beam under axial load 
are plotted in Figures 15 and 16. Similarly, deflections for transversal load case displacements at 
the free end of the beam are reported in Figures 17 and 18. 

Displacements from the interface element model compare very well to those from the 
conventional finite element model. Only in the transversal load case, a single interface element 
appears to slightly overestimate the vertical displacements. 

4. 4 Clamped- Clamped A symmetric Beam 

In this section an asymmetric beam, clamped at both ends, is loaded with two 
concentrated loads applied at two different points of the mesh. This problem was studied in 
[4,13] to validate the proposed interface method. The geometry of the problem and the position 
of the two concentrated loads are shown in Figure 19. 

The beam is meshed using two domains with different material properties. The ratio of 
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the Young’s modulus in the domains a and b is 0.2. The domain a is discretized with an 8 by 8 
mesh; two different meshes are adopted for the domain b. The mesh used for the conventional FE 
analysis is set to be 8 by 8, so that at the interface there is the same number of nodes. A different 
mesh is needed for discretizing domain b in order to test the interface model’s capabilities. An 8 
by 12 mesh has been chosen. The two forces, having the same magnitude F=1 but different 
directions, produce a complex state of deformation with strong stress concentrations at the 
interface that can challenge the interface element. 

Figures 20 and 21 compare the deformed configuration, produced by a conventional 
FEM, with that obtained using a FEM with four interface elements at the interface. No visible 
differences can be observed. 

The test is completed by a comparison between the in-plane displacements u and v along 
the interface as shown in Figures 22 and 23. The four interface elements connect two interfaces 
with a different number of nodes. Both displacement components are in satisfactory agreement 
with the conventional analysis. 

4.5 Simply Supported Plate Under Sinusoidal Load 

A simply supported plate is subjected to a transverse load distributed over the surface 
according to the following expression: 

7TX . Tty ~; n 

q = q Q sm — sin — 
a b 

In which q 0 represents the intensity of the load at the center of the plate, a and b are the plate 
dimensions, respectively, in the in-plane directions x and y. 

We assume the plate to be square with equal side dimensions, a = b = 2 . The value of q 0 
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is taken to be 1 . Taking advantage of geometrical and loading symmetries, only a quarter of the 
plate is analyzed. As shown in Figure 24, a fine mesh is used for a square near the center, while 
the remaining surface adopts a coarse mesh. The four noded Abaqus shell element S4R is 
employed to mesh both domains. Two interface elements on each edge are used to connect the 
two domains. Others properties are: E ] = E 2 = \e6 , v = 0.25 and thickness = 0.1. Figure 25 plots 
deformed and undeformed configurations of the plate as seen from the bottom. 

A classical plate solution to this problem exists in the literature. In Figure 26, the 
transverse displacements w along x at y— 1 are reported. Even in the presence of a high gradient 
of deformation, the model behaves well. This further confirms our confidence that the developed 
method is both robust and accurate. 
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5. CONCLUSIONS 


In the present work, an effective and robust interface element technology has been 
developed using the penalty method. This approach overcomes the numerical difficulties 
associated with the existing methods based on Lagrange multipliers. Additionally, the present 
formulation leads to a computational approach that is very efficient and completely compatible 
with existing commercial software. Significant effort has been directed toward identifying those 
model characteristics (element geometric properties, material properties and loads) that most 
strongly affect the required penalty parameter, and subsequently to developing simple 
“formulae” for automatically calculating the proper penalty parameter for each interface 
constraint. A wide variety of one and two-dimensional problems has been investigated 
analytically and computationally in order to correctly identify the analytical functional form of 
the penalty parameter for each constraint within many classes of problems. 

The resulting interface element is particularly efficient for finite element modeling of 
composite structures. When the material properties vary across and/or along the interface, the 
present method is often much more accurate than adopting a unique value of the penalty 
parameter. 

The present interface element has been implemented in the commercial finite element 
code ABAQUS as a User Element Subroutine (UEL), making it convenient to test the behavior 
and accuracy of the interface element for a wide range of problems. This approach has been 
validated by investigating a variety of two-dimensional problems, including composite 
laminates. 
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APPENDIX 


The types of finite elements that have been investigated are: conventionally formulated 
and reduced integrated Timoshenko beam elements, plane stress quadrilateral elements and plate 
elements based on the first order shear deformation theory (FSDT), or Mindlin plate theory. For 
each finite element formulation, different relations between penalty parameters and model 
characteristics (element geometric properties, material properties and loads) have been found. 
Results of this investigation are briefly presented in this appendix. 

Reduced integrated Timoshenko beam element 

Consider the problem of two linear reduced integrated Timoshenko beams connected 

through an interface element. The interface node is identified as i, the last node of the first beam 
is identified as 1 and the first node of the second beam is identified as 2. The Timoshenko beam 
element has three independent nodal DOFs: the axial displacement u, the transverse displacement 
w and the rotation y/. Thus, three different penalty parameters y n , y w and y v are employed to 

enforce the interface continuity constraints on the DOFs u, w and y/. In Table A1 are defined the 
penalty parameters to be associated with each of the constraints. 

Results obtained in our research lead us to affirm that the penalty parameters should be 
related to model characteristics as shown in Table A2. Where two values are present for the same 
DOF the bigger is chosen. Moreover, the subscript 1 and 2 marks, respectively, the material and 
geometrical properties belonging to the first and to the second beam. 

Plane stress quadrilateral elements 

As for the Timoshenko beam case, consider the problem of two plane stress quadrilateral 
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elements connected through an interface element. The generic node on the left of the interface is 
named 1, while the one on the right is named 2. Each node has two DOFs: the in-plane 
displacements u and v. Under conditions of axial load (load perpendicular to the interface 
surface), the following expression for the penalty parameter has been derived: 

v t \ 

n =P- — t=i.2 

V a k 

Where t fc and afc are, respectively, the thickness and the width of the element. 

This expression is very simple, but still contains information that would be difficult to 
obtain, namely: the width of the element, a. Inside the UEL subroutine only the coordinates of 
the nodes on the interface are known. Thus, assuming good modeling practice in FE of using 
elements with aspect ratio close to unity, we choose to approximate the dimension in the 
direction perpendicular to the interface with the one parallel to it. According to the adopted 
symbols, the expression for y k becomes: 



This simplification strongly affects the expression obtained for the flexural load case 
(load parallel to the interface surface), which is modified in the following way: 


Yk = 


' [ >;+ a 1 )g,/T 

4a, (2a, 2 i h ; ) _ 




A final step is required in order to obtain a single expression for an automatic evaluation 
of the penalty parameter. It consists in eliminating the scalar 4 from the denominator of the last 
expression for y k . This simplification is important since the interface could have any orientation, 


so it is better to not make any distinction between horizontal and vertical displacement DOFs. In 
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conclusion, only two penalty parameters are needed: y , for both the DOFs u and v of the nodes 
on the left of the interface, y 2 for the ones on the right of the interface. 


Y\ =P 


r Et ^ 
fri 

V b \ J 


Yi = P- 


( E t \ 
£ j 2 1 2 

K P ) 


Plate elements 

In plate elements each node has six DOFs: the in-plane displacements u and v, the 
transverse displacement w and the rotational DOFs (0 X , 9 y , 0,). The in-plane, transversal and 
rotational DOFs require different penalty parameters for optimal behavior. In particular: one, y u , 
is assigned to the in-plane DOFs (w, v), the second, y„ , to the transverse DOF w and the third, 
y g , is shared by all the rotational DOFs (0 X , 9 y , 0 Z ). The constraints are associated to the 

expression for computing the penalty parameters as in Table A3. 

Results obtained using the presented method have been simplified as for the plane stress 
quadrilateral element by approximating the element dimension in the direction perpendicular to 
the interface with the one parallel to it. The simplified relations are shown in Table A4. 
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T able 1. Tip axial displacement numerical results for beam extension. 



Mesh 

Solution 

Abaqus 

20x8 

4.000 

Abaqus with interface Element 

10x8 Left / 10x2 Right 

4.000 


10x8 Left / 10x4 Right 

4.000 

Classical 

- 

4.000 


Table 2. Tip deflection numerical results for beam flexure. 



Mesh 

Solution 

Abaqus 

20x2 

258.800 

20x4 

260.000 


20x8 

260.400 

Abaqus with interface Element 

10x8 Left / 10x2 Right 

259.500 

10x8 Left / 10x4 Right 10x8 

261.100 


1 0x8 Left / 1 0x8 Right 

260.400 

Classical 

- 

265.6 


Table Al. Penalty parameters associated with each of the constraints enforced in the Timoshenko 
beam element. 


Penalty parameter 

Constraint enforced 

Y„\ 


Y, i 

(w, -w,) 

Y v i 


Yu 7 


Y w2 

(w,-w 2 ) 

Yv 2 

(V.-Vl) 
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Table A2. 









Table A3. 


Penalty parameters associated with each of the constraints enforced in a Mindlin plate 
element. 


Penalty parameter 

Constraint enforced 

Y„ i 

(«, - »i ) 

r,i 

( v / - v . ) 

r w i 

(w,-w,) 



Yu 2 

(u,-u 2 ) 

fv2 

( v /- v 2 ) 

Y„2 

(w,-w 2 ) 

Y qi 

(<?-«.) 
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(b) 

Figure 2. (a) Beam element connected to a fixed point V by one 
interface element; (b) two linear beam elements connected in 
Y using an interface element and clamped at one tip. 
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Figure 5. Stability of the solution for different values of y- bending 
load. 
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Figure 8. 


(b) 

(a) Interface element and (b) conventional FE mesh - zoom of 
the deformed configuration. 
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S22 VALUE 



-1 . 01DE+00 

- 8 . 9 2 4 E - 01 

- 7 . 7 5 2 E - 01 
-6 . 580E-01 
-5 . 407E-01 
-4 .235 E-01 
-3.063E-D1 
-1 . 89 0 E-01 
-7.181 E-02 
+ 4 . 542E-02 
+ 1 . 62 6 E-01 
+ 2 . 799 E-01 
+3 .971 E-01 
+5.143E-01 



S22 VALUE 



-1 .01 3E+00 
-8 .950E-D1 
-7 . 7 7 5 E-01 
-6 . 6 0 1 E-01 
-5 . 4 26 E-01 
-4 . 2 51 E-01 
-3 . 076E-01 
-1 .901 E-01 
-7 . 263E-02 
+ 4 . 4 86 E-02 
+1 . 62 3 E-01 
+ 2 .798 E-01 
+ 3 . 973 E-01.. 
+ 5 . 1 4 8 E-01 



(b) 


Figure 10. (a) Interface element and (b) conventional FE solution - 
stress distribution in the direction 2 (a 22 ). 
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U1 VALUE 

- +6 . 956E-38 

- +7.717 E+00. 
- +1 .543E+01 

+2.31 5E + 01 
+3.0B7E+01 
+ 3.859 E + Dl 
+ 4.630 E + 01 
+5 . 402E+01 
+6 . 174E+01 
+ 6.946 E + 01 
+7.717 E+01. 
+ 8 . 489E + 01 
+9.261 E+01 

- +1 .003E + 02 





VALUE 
+1 . 54BE-37 

+7 .71 7E+00 
+1 . 543E+01 
+2.31 5E+01 
+3 . 0 8 7 E + 01 
+ 3 . B 5 9 E + 01 
+4 . 630E+01 
+5 .402 E + 01 
+ 6.174 E+01 
+ 6.94 6 E+01 
+ 7.717 E + 01 
+ 8 . 4 89 E + 01 
+9.261 E+01 
+1 .003E+02 



Figure 11 . (a) Interface element and (b) conventional FE solution - 
horizontal displacement distribution (Uj. 
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VALUE 

■ -1 . 5 3 3 E + 01 

-- 1 .415E+01 

- -1 . 297E + 01 
-1 .179E + 01 
-1 .061 E+01 
-9 .432E+00 
-8.253 E + OD 
-7 . 074E+00 
-5 .895E+00 
-4.716 E + 00 
-3.537 E+00 
-2 . 358E+00 

• -1 . 1 7 9 E+OO 

- +4 . 160E-38 



VALUE 

• -1 . 533E + 01 
■ -1 . 4 1 5 E+01 

-1 . 297E+01 
-1 . 179E + 01 
-1 .061 E+01 
-9 . 4 33E + 00 
-8 . 25 4 E + 00 
-7.07 4 E+00 
-5 . 89 5 E + 00 
-4 .716 E + DO 
-3 . 537E+00 
-2 . 3 58 E + 00 
-1 .179 E + 00 

• +3 .821E-38 



I 


Figure 12. (a) Interface element and (b) conventional FE solution 




vertical displacement distribution (U 2 ). 






Abaqus Conventional 



Figure 15 . Beam extension along the thickness at the free end under axial load. 




Abaqus Conventional 



Figure 16 . Beam vertical deflection along the thickness at the free end under axial load. 



Abaqus Conventional 



Figure 17 . Beam extension along the thickness at the free end under transversal load. 




Abaqus Conventional 



Figure 18. Beam vertical deflection along the thickness at the free end under transversal load. 






Os 
» -n 


w 



c 

o 


o 

m 


w 

(Jh 

in 

3 

a" 

a3 


C 3 

e 

o 


s 

§ 


o 

o 



*T 3 

0 ) 

S3 

*3 

o 

S3 

O 





O 

o 


T 3 

<U 


£ 

<L> 


O 

<N 


4 > 

L< 


s 

ftD 

to 


o 

VO 


JC 

vi 



Figure 21 . Deformed configuration obtained using interface elements approach. 


Abaqus Conventional 



Figure 22. Axial displacements along the thickness at the interface. 



Abaqus Conventional 



Figure 23. Transversal displacements along the thickness at the interface. 
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Figure 26. Variation of the transverse displacement w along x at y = 




